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I. INTRODUCTION 


Q 


This paper is devoted to a study of the nonautonomous Adler equation 

9 = r(t) — sin 9. (1) 

When r is independent of time this equation describes phase synchronization between a pair 
of coupled oscillators. In this case 9 = (j)i — (j )2 represents the difference in the phases (j)j 
of the two oscillators and r represents the normalized frequency difference. When |r| < 1 
the equation describes a phase-locked state; when |r| > 1 the phase difference increases 
or decreases monotonically, corresponding to repeated phase slips. The transition between 
these two states is an example of a SNIPER (saddle-node infinite period) or SNIC (saddle- 
node on an invariant circle) bifurcation In this bifurcation the phase slip period diverges 
like \l\/r — \ as r decreases towards r = 1, m contrast to transitions associated with global 
bifurcations. 

The nonautomous equation ([1]) with r = r{t) and r{t) a periodic function of time thus 
describes the effects of temporal modulation of the SNIPER bifurcation. Such modulation is 
of interest since for part of the modulation cycle the oscillators may be phase-locked while for 
the rest of the cycle they may undergo phase slips. In this paper we show that the interplay 
between these two states is complex, and characterize the resulting behavior for both high 
and low frequency modulation r(t); the intermediate case in which the modulation period 
is comparable to the phase slip period is of particular interest and is also investigated here 
in detail. 

The nonautonomous Adler equation arises in a number of applications. First and foremost 
it arises in systems of driven identical active rotators js, or, equivalently, driven arrays 


of Josephson junctions 




described by the equations 

M 


(j)j = oj(t) — sin (f)j — K sin(0j 


( 2 ) 


rn=l 


Here u is the intrinsic frequency and K measures the coupling strength. In terms of the 
Kuramoto order parameter, i?expi$ = this system can be written in the 

equivalent form 

9j = uj{t) — a — KRsm9j, (3) 

where 9j = (pj — a, KR = a/1 + {KRy + 2KRcosQ and tana = KRsin ^{1 + KRcos <f))“^. 
Since R and $ are in general functions of time y] the quantities R and a will also be 
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functions of time and these are determined by the collective dynamics of the remaining 
M — 1 oscillators. When M is large the latter are unaffected by the behavior of an individual 
oscillator, and R and a can therefore be assumed to be given. The dynamics of each oscillator 
are thus described by an Adler equation with a time-dependent effective frequency and a 
time-dependent effective coupling constant. The latter dependence can be removed using 
the simple substitution dr = KRdt provided K{t) remains bounded away from zero. 

Nonautonomous effects also arise in phase-coupled oscillator systems of Kuramoto type 
and these are of interest in neural models. In models of this type the coupling strength 
Kjk between oscillators j and k is taken to be a function of time, reflecting either evolution 
12 1 or the effects of a drug, during anesthesia, for example [l^. The 


of the network 


simplest model of this type. 


M 


= Ul 




( 4 ) 


m=l 


can be written in the equivalent form 


9j = u — ^ — K{t)R{t) sin 9j, 


( 5 ) 


where 9j = — <h. Thus the dynamics of each individual oscillator are determined by the 

global behavior of the system through the quantities KR and $. When M is large both R 
and $ may be taken as given, independent of the behavior of the oscillator j. The resulting 
system can be cast in the form 

9j = u{t) - sm9j, ( 6 ) 

where the prime denotes differentiation with respect to r, dr = KRdt and a)(r) = 
[uj/K{t)R{t)] — with K{t) = K[t{T)], R{t) = i?[t(r)] etc. It suffices, therefore, 

to consider the effects of a time-dependent effective frequency only. Related models arise in 


systems with frequency adaptation 


14i\ . An excellent review of the origin of nonautonomous 




effects in the Kuramoto model and its variants can be found in |l5l |. 

Finally, the nonautonomous Adler equation also describes a single resistively shunted 


Josephson junction driven by a biased AC current jl^ . 
equation, motivated by observations of Shapiro steps jlT 




heoretical investigations of this 
in the supercurrent, have illu- 


mmated a wealth of mode-loekiag behavior MM- Large arrays of coupled Jos^^on 

junctions are thus amenable to the same type of analysis as active rotator systems [a, l2l[ |. 
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The paper is organized as follows. In the next section we snmmarize the basic properties 
of the Adler eqnation with and withont time-periodic modnlation. In Sec. Illll we stndy, nnder 
a variety of conditions, periodic orbits of the nonantonomons Adler eqnation that take the 
form of oscillations abont a phase-locked state. In Sec. HV] we stndy the so-called phase¬ 
winding trajectories describing repeated phase slips and identify the regions in parameter 
space where different states of this type are found. In Sec. |V] we show that an adiabatic 
analysis describes accurately the resulting parameter space not only for low modulation 
frequencies but in fact remains accurate far outside of this regime. Section |Vl] provides a 
brief summary of the results and discusses potential applications of the theory. 

II. THE ADLER EQUATION 

The Adler equation ([T]) with constant r has several symmetries of interest. The equation is 
invariant under complete rotations W : 6* —?■ 6+271, and time translations 7^ : t —)■ t-|-r by an 
arbitrary real r. In addition, it is invariant under the phase symmetry Vq : {t, 9) —)■ (—t, tt—O) 
and the parameter symmetry TZq : {r,d) —{r,6). As already mentioned, the fixed points 

or equilibria of Eq. ([1]) correspond to phase-locking between the two oscillators, and these 
exist in the parameter interval |r| < 1: 


9eg = sin ^ r. (7) 

If 6 is defined mod 27r, this condition determines two branches of equilibria that merge 
in a saddle-node bifurcation at r = ±1 and are related by Vo- One of these branches is 
stable and can be identified by the condition dr6eq > 0 while the other is unstable and is 
characterized hj drdeg < 0. No fixed points exist for |r| > 1: 6 increases monotonically when 
r > 1 and decreases monotonically when r < — 1. When 6 is defined mod 27r the resulting 
trajectories are both periodic in time and the steady state SNIPER bifurcations at r = ±1 
generate periodic orbits, a consequence of the global organization of the stable and unstable 
manifolds of the hxed points. 

In the present work we find it convenient to think of 0 as a variable defined on the real 
line. When this is done the equation has an inhnite number of stable and unstable equilibria 
that differ in the number of 27r turns relative to an arbitrary origin 6 = 0. We refer to these 
turns as phase slips since one of the two oscillators is now ahead of the other by an integer 
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number of 27r rotations. Trajectories outside of the phase-locked region will incur positive 
or negative phase slips with frequency 

ojq = \/r^ — 1 . ( 8 ) 

This frequency approaches zero in a characteristic square root manner as |r| approaches 
|r| = 1 from above j^. 

When the frequency parameter r oscillates in time, 


r = ro -I- asin(27rf/T), 


(9) 


the system retains the winding symmetry W, while the translation symmetry becomes 
discrete T : f —)■ f -f T. The phase symmetry now includes a time shift, V : (f, 9) —)■ 
(T/2 — t,7i — 9). The parameter symmetry takes the form 71 : (ro, a, 9) —(ro, a, 9). There 

is also an additional parameter symmetry S : {t,a) (t -|- T/2, —a). We remark that, as 
already explained, any time-dependence in the coupling parameter K > 0 can be removed 
by a simple transformation, and this parameter is therefore scaled to unity. 

Depending on the amplitude a and the period T of the frequency modulation ([9]) the 
solutions of the resulting equation take the form of oscillations about a phase-locked state 
or describe repeated phase slips in which the phase difference 9 drifts with a nonzero mean 
speed. We identify below a series of resonances between the modulation period T and the 
time scale for the generation of a phase slip. The resulting parameter space structure is 
determined using a combination of numerical simulations, numerical continuation 22| and 
asymptotic methods. Regions with an integer number of phase slips per period are separated 
by regions with noninteger numbers of phase slips, and include canard trajectories that drift 
along unstable equilibria. Both high and low frequency modulation is considered. We do 
not consider noise-triggered phase slips. 


III. PERIODIC ORBITS 

Phase-locked states of the autonomous system ([T]) may undergo phase slips in the pres¬ 
ence of modulated frequency while remaining phase-locked on average. For such solutions 
the number of negative phase slips balances the number of positive phase slips over one 
modulation cycle. Figure [1] shows the bifurcation diagram for the nonautonomous Adler 
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(a) (b) 




^0 r 

FIG. 1. (Color online) (a) Bifurcation diagram showing the average phase (9) = T~^ Jq 9{t)dt 
of periodic orbits as a function of vq when a = 2 and T = 15 (blue dashed line), T « 23.01 (red 
dash-dotted line) and T = 25 (black solid line), (b) Sample trajectories, in corresponding line type, 
in the (r, 6) plane for solutions with rg = 0 and (6) = 2tt and 7tt, superposed on the branch of 
equilibria of the autonomous system (a = 0), represented by a green dotted line. 

equation ([1]) with the periodic modulation ([9]) along with sample trajectories at two points 
on the solution branches, both superposed on the corresponding equilibrium solutions of the 
autonomous system, i.e., r = rg. The solution branches snake, i.e., they undergo repeated 
back-and-forth oscillations as the parameter rg varies. The extrema of these oscillations 
correspond to the SNIPER bifurcations at r = ±1; the equilibria with a positive slope cor¬ 
respond to stable solutions while those with a negative slope are unstable. Thus along the 
branch of equilibria stability changes at every fold. 

The trajectories shown in Fig. [I](b) are periodic, with period T, and their bifurcation 
structure parallels that of the phase-locked states in the autonomous system: the solutions 
snake within an rg interval determined by a pair of folds on either side as shown in Fig. [I](a). 
The amplitude of this oscillation and its shape depends on the period T of the forcing which 
also affects the solution stability. For example, for (6) = 2ti and rg = 0, the solution of the 
autonomous problem is stable, but becomes unstable for T = 15 as most of the periodic 
orbit tracks the unstable branch of the autonomous problem, before becoming stable again 
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FIG. 2. (Color online) The fonr distinct orbits generated on applying the symmetries {I,TZ,V,S) 
to the stable periodic orbit compnted for T = 15, tq = 0.2, and a = 2. A sequence of orbits with 
9^9 + 27rn can be found by applying W"' to each of the four solutions. These orbits lie on the 
branch displayed in Fig. UKa) for T = 15. The symmetry W has been applied in order to prevent 
overlap between the four distinct orbits. The equilibria of the autonomous system (a = 0) are 
shown as a green dotted line. 


for T = 25. A numerical computation of the Floquet multiplier exp y— cos 9{t)dij for 
the Adler equation linearized about the periodic orbit during the continuation procedure 
conhrms that the upward (downward) sloping portions of the solution branch remain stable 
(unstable) all the way to the folds. 

The presence of the symmetries allows to generate other solutions from the one calculated. 
Figure |2] shows the four different orbits produced by applying the eight different symmetries 
generated by (X, 77, P, 5): X, 77, P, 5,77P, 775, PiS, 77PiS to a periodic orbit obtained for 
tq = 0.2, T = 15 and a = 2. These periodic orbits he on the same solution branch in 
Fig. UKa). The symmetry S acts like the identity, the time shift compensating for the 
reversal of a. Application of T does not produce new orbits, and we can shift any periodic 
orbit to higher or lower values of 6 by multiples of 27r using powers of W. We take advantage 
of the latter to avoid overlap among the different solutions shown in Fig. |2l 

Figure [3] shows how the existence region of the periodic orbit, labeled PO, evolves with T. 












FIG. 3. (a) Locus of the folds that define the boundary of the PO region in the {ro,T) plane. 
The horizontal dashed and solid lines indicate the values of T corresponding to the branches of 
periodic orbits computed in Fig. HJ^a). (b) The amplitude A0 = ^max — ^min of a periodic orbit 
with ro = 0 and o = 2 as function of the period T. The dotted horizontal lines correspond to 
the pinched zones at T 9.33, 23.01 and 37.31 in panel (a); at these the corresponding periodic 
orbits are characterized by A0 ~ 4.95, 11.32 and 17.71 and deviate from multiples of 27r by 
(27rn — A0)/27r ~ 0.21, 0.20, 0.18, respectively. 


Numerical continuation of the folds at the two edges of PO reveals a series of pinched zones 
in which the folds “cross” and the PO region is reduced to the single value tq = 0. This 
accounts for the switch in the orientation of the branch as T increases (see Fig. [11(a)). We 
call the regions between the pinched zones sweet spots. Within each of these sweet spots, the 
number of positive and negative phase slips during one cycle is the same, and the orbits are 
therefore qualitatively similar. The resulting structure, shown in Fig. [3](a), is reminiscent of 


the structure observed in 


23| . Figure |3](b) shows the amplitude of the oscillation in 6 for 


periodic orbits at ro = 0 as a function of the period T. The hgure reveals that N positive and 
negative phase slips can occur even when A9 = 6*max — < ^nN. This is a consequence of 

the fact that the two successive saddle-nodes at r = ±1 are separated by a phase difference 
of TT. Figure m shows a series of periodic orbits that transition from zero to one positive and 
one negative phase slip as a (equivalently T) increases. 
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FIG. 4. (Color online) A series of periodic orbits (solid black) for T = 25, tq = 0 and increasing 
values of a, corresponding to increasing oscillation amplitude A0 = vr, Svr/d, 37r/2, 77r/4, 27r, super¬ 
posed on top of the bifurcation diagram of the phase-locked solutions of the autonomous system 
a = 0 (green dotted line). The transition from zero phase slips to one positive and one negative 
phase slip is indicated by a dashed blue line and corresponds to a « 1.29 and A6I ps l.bSvr. 

A. Birth of periodic orbits 

To understand the effect of a time-dependent frequency parameter on the dynamics of 
phase-locked oscillators, we start out by considering the high-frequency modulation limit of 
the Adler equation ([T]) with the time-periodic modulation ([9]). We write T = 27re/a;, where 
e <C 1, and u ~ 0(1) is a suitably scaled frequency, and dehne the fast time 0 by cut = e0. 
The Adler equation becomes 

udff)9 = e(rQ + asm.(l) — smO — dt9). ( 10 ) 

We assume that 0(0, t) = 0o(05 t) + e^02(0, t) + ... and carry out the calculation 

order by order. The leading order equation shows that 0o = 00 (^) is independent of the 
fast oscillation time. The 0(e) equation yields, after integration over the fast period of the 
forcing, 


5i0o = - sin 00 - 


( 11 ) 
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Thus, at leading order, the averaged system follows an autonomous Adler equation with 
constant forcing equal to the average of the periodically modulated case. 

The solution at order e reads 

9i{(j),t) = -—cos (j) +'ipiit), (12) 

OJ 

where ipi is determined through the solvability condition at the next order. This next order 
equation reads 

ujd ^62 = -9i cos 00 - dtOi, (13) 

and integration over the fast period gives the solvability condition 

dtipi =-i)i cos iiQ. (14) 

The solution at order is thus 

6 * 2 ( 0 , t) = sin (j) cos + 02 ( 6 ), (15) 

while the order equation reads 

= -62 cos 00 + | 0 ? sin 0 o - 8162, (16) 

leading to a solvability condition for 02 ^ 

0*02 + 02 COS 00 = ^ sin 00 + 101 sin 0o. (17) 

To study the average dynamics, we dehne the period-averaged phase 

p 27 V 

'll) = (27r) ^ / (00 + 601 + 6 ^ 02 ) d(j). (18) 

Jo 

This expression is accurate to order 0{e‘^). Summing the solvability conditions now yields 
the equation 

0*0 = ro - (1 - fS) sin0 + 0{T^), (19) 

where we have replaced cn/e by 27r/T. Thus, in the high-frequency limit, the averaged 
dynamics follows an Adler equation for which the amplitude of the nonlinear term that 
characterizes the coupling strength between the two oscillators decreases in proportion to 
(aT)^. The phase-locked region of the averaged equation (IT^ that defines the PO region for 
the time-dependent Adler equation thus exists for |ro| = 1 — (aT/dvr)^, and the introduction 
of high-frequency modulation narrows the width of the phase-locked region in the parameter 
ro by 2{aT/4'irf. 
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B. Death of periodic orbits 


Asymptotic analysis near the folds that dehne the edges of PO can provide some insight 
into the break-np of the periodic orbits. We consider pertnrbations abont the marginally 
stable orbit at the left (ro = r_) and right (ro = r+) edges of PO for a given modulation 
frequency to = 27i/T and amplitude a, namely Eq. ([T]) with r = ro + a sin cat, where ro = 
r± + and e <C 1. We use multiple time scales by introducing a slow time r = et on which 
the system incurs net phase slips and expand the phase variable as 9 = 9o + e9i + e‘^92 + ■ ■ ■ ■ 
The leading order equation, dt9Q = rj- + asincat —sin 6 'o, is solved by the marginally stable 
periodic orbit, which we have computed numerically via continuation. The 0(e) equation is 


dt9i + 9i cos 9o = —dr9o 


( 20 ) 


which has a solution of the form 9i = A exp (— J cos9Qdt^ for a slowly-varying amplitude A 
as 6*0 does not depend on the slow time. At the equation reads 


9 j 6*2 -I- 6*2 cos 9o = fi+ ^ 6*1 sin 9o — dr9i. 


( 21 ) 


The existence of a solution in 92 that is T-periodic requires that the solvability condition 


dr A — fiai -\- — 0 ( 2 ^^ 


( 22 ) 


be satisfied where the coefficients can be computed numerically from the integrals 


T 


ai = — exp ( / cos 6 *o(it ] dt, a 2 = — sin 6 *oexp ( — / cos d^dt ) dt. (23) 


T 


Thus, just outside of PO, the system will incur net phase slips with a frequency 


^^Eiip ^ ,/2 |qiQ2(j’o - J’±)|- 


(24) 


Figure O shows a comparison of this frequency as a function of tq with simulations near the 
right edge of PO for T = 15, where r_|_ 0.305, and a = ^/2\oo(^\ k, 1.163. The coefficient 

that describes the square root dependence of the frequency on the distance from the left 
edge of PO will be identical to the one computed for the right edge owing to the symmetry 

7 ^. 
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FIG. 5. (Color online) (a) A plot of the freqnency flgiip at which phase slips occnr just outside of 
PO as a function of the distance from the edge when a = 2 and T = 15. The solid green 

line is the prediction in Eq. (I24|] from asymptotic theory while the blue dots are computed from 
time simulations. 


C. The asymptotics of sweet spots 

When large excursions of the forcing parameter are allowed during a high-frequency cycle, 
a balance is struck that allows a finite number of phase slips to occur. We keep T = 27re/u 
but link the amplitude of the forcing to the frequency hy a = p/e = 2np/uT. Upon dehning 
the fast time-scale ej) = u}t/e, the Adler equation becomes 

ojd^O — psmcf) = e{rQ — smO — dtO). (25) 

Using an asymptotic series of the form = 00(0, f) -|- e0i(0, f) -f e^02(0, t) + • • • and 

solving the leading order equation we obtain 

eo{(j),t) =-^cos(p + ^/Jo{t). (26) 

The evolution of 00 is determined from a solvability condition at next order. Since the order 
e equation reads 

ud^di = ro -b sin cos 0 - 0o) - dti/o, (27) 

the required solvability condition is 


dt<l’a = ro- Jo(PsmV>o, 


(28) 
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(a) (b) 




FIG. 6. (a) The PO region in the {ro,aT) parameter plane corresponding to stable phase-locked 

solutions of the Adler equation when the forcing has high frequency and a large amplitude, (b) 
The leading order amplitude A9 = Omax — dmin of a periodic orbit at tq = 0 as a function of aT /27r. 
Horizontal dotted lines correspond to the first three pinched zones which coincide with the zeros 
of Jq: aTj2TT « 2.40, 5.52 and 8.65. 

where Jq is the Bessel function of the first kind. The averaged dynamics thus follow 
an autonomous Adler equation with a constant frequency and a coupling strength given 
by Jo{p/u) = Jo{aT/27i). The boundaries of the PO region are thus dehned by tq = 
±Jo(aT/27r) and these oscillate in ro as aT increases with an amplitude that decreases with 
increasing T (Fig. |6]). The location of the pinched zones is thus determined by the zeros 
of jQ{aT/27i). Between these are the sweet spots where periodic orbits exist over the hnite 
range |ro| < \Jo{aT/27r)\. The reversal of orientation of the folds seen in Fig. [T](a) is anal¬ 
ogous to sign changes of Jo{aT/27r) in this high frequency, large amplitude limit, as shown 
in Fig. El 


D. Amplitude dependence 

We now examine how periodic solutions within PO behave as a function of the amplitude 
of the modulation by hxing tq = 0 and performing numerical continuation in a (Fig. [7]). As 
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FIG. 7. Bifurcation diagrams showing (a,c) the average phase (0) = T~^ 6{t)dt (solid lines) and 
(b,d) the oscillation amplitude A0 = 0max — ^min of periodic orbits as a function of a when ro = 0 
and T = 25. The solutions shown in (a) collapse onto a single curve when plotted in terms of A0 
in (b). When ro = 0.1 and T = 25, the grid structure of (a) separates into isolated loops shown in 
(b) that collapse onto disconnected line segments when plotted in terms of A0 in (d). 


long as ro is in the interior of PO, each value of a admits two periodic orbits on a 27r interval 
for {6). One is stable, one is unstable, and they are related by the phase symmetry V. The 
symmetries of the system further imply that the locations of these orbits at ro = 0 are hxed 
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FIG. 8. (A)-(C) Periodic orbits with {6) = tt in the (r, 6) plane when rg = 0, T = 30 and 

a = 1, 1.5, 2. (C)-(E) Periodic orbits with {9) = it, 27r, Svr in the {r,9) plane when ro = 0, T = 30 
and a = 2. The orbits correspond to the red dots in Fig. ma) labeled with capital letters. 

at (6) = rmi for m G Z and such solutions persist for all values of a (horizontal lines in panel 
(a) of Fig. [7]). The pinched zones where the PO boundaries cross (Fig.[3](a)) and the snaking 
branch becomes vertical (red dash-dotted line in Fig. [T](a)) correspond to codimension two 
points in the (a, T) plane; at these points a continuum of periodic orbits parametrized by 
the phase average (6) is present. Thus when tq = 0 the periodic orbits create the grid-like 
bifurcation diagram shown in Fig. [7](a). This grid structure breaks apart into isolated loops 
of solutions as soon as ro ^ 0 , and gaps between the regions of existence of periodic orbits 
begin to emerge (cf. Fig. [3](a)). The loops that emerge from the breakup of the rectangular 
grid structure at ro = 0 when ro 7 ^ 0 shrink to zero with increasing a (or T), as expected 
from Fig. [3](a). Numerical continuation of the boundary of the PO region as a function of 
a when vq = 0.1 and T = 25 reveals that periodic orbits persist only to a ~ 14.5. 

Figure [8] shows solutions for ro = 0 with parameters values indicated in Fig. [7](a) by red 
dots labeled with the corresponding capital letter. The equilibria for the autonomous prob¬ 
lem are shown for reference (dotted line). These reveal that the periodic orbits alternately 
track branches of unstable and stable equilibria for part of each oscillation cycle (orbits A, 
B, C), and likewise for C, D, E. Since orbits that track stable equilibria are expected to be 
stable when the drift along such equilibria is sufficiently slow, we expect that orbits B and D 
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FIG. 9. Periodic orbits along the first vertical solution branch in Fig. [7] in the (r, 0) plane when 
ro = 0, T = 25 and a ~ 1.2. These solutions are characterized by {9) that is a fraction of 27r, viz. 
TT, 57r/4, 37r/2, Tvr/d and 27r. The orbits correspond to the unlabeled blue dots in Fig. [7[a). 

are stable while A, C and E are unstable. This expectation is confirmed by explicit stability 
calculations. 


E. Canards 

Figure [9] shows periodic orbits from the hrst vertical solution branch in Fig. [TJ^a) cor¬ 
responding to the dark blue dots not labeled with capital letters. These periodic orbits 
all have the same value of a ~ 1.2 and correspond to pinched zone solutions with {9) = vr, 
57r/4, ?)'k/2, 77r/4 and 27r. These solutions illustrate how the periodic orbit expands to larger 
{9) while tracking the equilibria of the autonomous system. These are beginning to reveal 
characteristics of the so-called canard states familiar from studies of slow-fast systems. For 
example, the third panel shows an orbit that slowly tracks a branch of stable equilibria to¬ 
wards lower 9 and smaller r followed by tracking a branch of unstable equilibria towards yet 
smaller 9 but increasing r, before an abrupt transition near the right fold that restores the 
original 9 value essentially instantaneously, i.e., at essentially constant r. This difference in 
timescales is not as pronounced in the last panel of Fig. |9]but can be enhanced by increasing 
the modulation period. Figure ITO] shows typical canard trajectories with a clear separation 
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(a) (b) 





FIG. 10. (Color online) (a) Two-headed canard trajectories 9{r) for tq = 0, T = 100 and 
a PS 1.064807, 1.064865, 1.064871, 1.064872, 1.064876, 1.066086, 1.177531 and 1.198182. (b) The 
corresponding solutions 9{t) and 9{t). 


of timescales, obtained for T = 100, tq = 0 and slightly different modulation amplitudes a. 
Increasing the amplitude a very slightly leads the canard to overshoot the right saddle-node 
and can make it depart from the branch of unstable equilibria upwards, i.e., in the opposite 
direction as compared to the solutions for slightly smaller a. The latter case leads to a 
different type of canard: the system jumps from the unstable solution branch to the upper 
branch of stable equilibria, which it then follows downward in 9. After reaching the upper 
left fold of the equilibria the trajectory jumps to the lower left fold and thereafter follows 
the lower unstable equilibria towards larger r, resulting in the same sequence of transitions 
but now as r increases. The resulting solution is periodic but is characterized by phase 
slips that take place inside the snaking region |r| < 1. This behavior is exemplihed by the 
outer canard trajectory in Figs. [TUf a): the associated 6 displays an inverse peak, as shown 
in Fig. [TUlf b). Additional time separation along the stable and unstable manifolds can be 
achieved by increasing T further. For example, Fig. [TT] shows several “two-headed” canard 
trajectories obtained for T = 300. 
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(a) (b) 





FIG. 11. (Color online) Two-headed canard trajectories computed by numerical continuation of 
periodic orbits in the parameter a. The parameters are ro = 0, T = 300 and a ~ 1.02115308, 
1.02116560, 1.02116562. 

IV. WINDING TRAJECTORIES 


Outside of the phase-locked region of the Adler equation with constant frequency param¬ 
eter (ro G [—1,1], a = 0) there exist winding solutions that complete phase slips with the 
frequency given by Eq. ([H]). The introduction of a modulation in r with period T (Eq. |9l 
a 7 ^ 0) generates winding solutions even when the average value tq lies within [—1,1]. This 
occurs for values of ro outside of PO (but |ro| < 1), and is a consequence of an imbalance 
between positive and negative phase slips. 

We dehne the winding number of a trajectory in the modulated system as the average 
number of net phase slips per period, 


N = lim 

m^oo 


e{mT) - e{o) 

27rm 


(29) 


with m G Z. Figure [12] shows solution branches with integer winding numbers N = 1,2,3 
when a = 2 and T = 25 (solid lines). These were computed by numerical continuation as 
a boundary value problem with the constraint that 0{T) — 6*(0) = 27rN. Trajectories with 
integer winding number exist over hnite ranges of the parameter ro. Solutions displaying 
an extra positive phase slip over each modulation cycle have winding number V = 1; these 
exist for ri min ~ 0.1 < ro < ri ~ 0.4. To the right of this interval he solutions with 
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(a) (b) 




FIG. 12. (Color online) (a) The phase (9) = T~^ 9(t) dt averaged over T of winding orbits as a 

function of tq when a = 2 and T = 25. Since 9 is no longer periodic all points with the same (9), 
mod 27r, at a particular value of tq lie on the same trajectory. The black (with circle), red (with 
square) and blue (with triangle) branches have winding numbers iV = 1, 2, 3, respectively. The 
branches of solutions with the same winding numbers but constant frequency parameter r = tq 
are shown as (vertical) dotted lines, (b) Sample winding trajectories corresponding to the colored 
symbols in panel (a). 

winding nnniber N = 2, extending from r 2 ,min ~ 0.4 to r 2 ,max ~ 0.6. Solutions with higher 
integer winding number exist beyond this point as exemplihed by the N = 3 solutions in 

Fig. [121 


A. Resonance tongues 


The parameter range containing integer winding solutions forms through the opening of 
resonance tongues as the modulation amplitude a increases from zero. We write, following 
ji^ . X = tan6*/2 to put the Adler equation in the form 


i. = -r — X + -rx . 


(30) 
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The Riccati transformation x = —‘lyjry now generates a second order linear equation for 
the variable y{t\- 

• \ 2 

+ = ( 31 ) 

_ 1 . r _ —dt 

Using the standard transformation y = ze ’■ we hnally obtain the Hill equation 






(32) 


Substituting the time-dependent frequency parameter r specihed in Eq. 
a -C 1 yields the Mathieu equation 


and assuming 


'~4~ ^ ^ + (’'o - sm{uit - 0 ) z + = 0. 


(33) 


where u = 27i/T and tan^ = (rg — u‘^)/u. Phase slips in the original Adler equation 
correspond to divergences of x; these in turn correspond to zero crossings of y and z. 

The resonance tongues grow in this asymptotic limit according to the characteristic curves 
of the above Mathieu equation. We compare these asymptotic predictions with the numerical 
computation of the resonance tongues through two-parameter continuation of folds on the 
branches of winding solutions. The tongues associated with the 1:1, 2:1, and 3:1 resonances 
between the winding frequency and the modulation frequency are shown in Fig. [13] alongside 
the predictions from the characteristic curves of the Mathieu equation (i33|) . 

The resonance tongues enter farther into the phase-locked region |ro| < 1 as a increases. 
We observe that as a increases the location of the tongues begins to depart from the Mathieu 
equation predictions, as noted already in the context of Josephson junction models [l6[ (Ch. 
11). In particular, the interaction of these tongues with their negative winding counterparts 
leads to qualitative changes: for a > 1.29, the width of the 1:1 resonance tongue stops 
growing monotonically and its left boundary turns abruptly from rg 0 to larger rg; at 
rg 0.25 , a ~ 1.57 the tongue collapses to a single point before growing again. This 
situation repeats as a increases and the tongue therefore describes a succession of sweet 
spots and pinched zones. The same behavior is observed for the subsequent resonance 
tongues: the 2:1 resonance tongue starts to shrink at rg ^ 0.25, a ~ 1.57 and collapses to a 
point at rg 0.50, a ~ 1.86, etc. 
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(a) (b) 



FIG. 13. (a) Resonance tongues for the 1:1, 2:1 and 3:1 resonances between the winding frequency 
and the modulation frequency in the (rg, a) plane when T = 25. The resonance tongues correspond 
to the solution branches shown in Fig. [T2]with 1, 2 and 3 phase slips per period of the modulation 
cycle, respectively. The boxed region in the lower right of panel (a) is replotted in panel (b) along 
with the predictions for the location of the tongues from Eq. (1331) in dashed lines. 

B. Partitioning of the parameter space 

The parameter plane (ro,T) can be partitioned in terms of winding nnmber by following 
the folds of the A^:l resonant winding trajectories snch as those shown in Fig. [121 The 
resnlting partitioning of parameter space is shown in Fig. [141 To obtain this fignre the 
branches with winding nnmbers 1 < iV < 7 were continned in rg for a = 2 and T = 5, 
followed by continnation of the saddle-nodes on these branches in the parameters rg and 
T. The region PO of periodic orbits was compnted in a similar way for completeness. The 
sweet spot and pinching strncture of regions with constant integer winding nnmber that 
begins to emerge in Fig. IT^ a) can also be seen as T increases for hxed a. The width of these 
sweet spots decreases with T. For inhnite periods, any small departnre from the symmetry 
axis rg = 0 leads to the dominance of positive or negative phase slips over the other. 

Tims the parameter plane is partitioned into regions with solntions displaying zero, one, 
two or more net phase slips per cycle. Each of these regions possesses a strnctnre similar 
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FIG. 14. (Color online) Average winding number per period T of the frequency parameter shown in 
the {rQ,T) plane for a = 2. No net phase slips occur over the course of a modulation period in the 
dark region to the left; the alternating lighter yellow/darker orange regions to the right indicate 
1,2,3,... net phase slips as tq increases. The (lightest) gray transition zones have non-integer 
winding numbers. Trajectories with negative winding number are located in regions obtained by 
reflection in rg = 0. 

to that of the PO region with zero net phase slips. The hrst region to the right of PO 
corresponds to solntions that nndergo one extra positive phase slip within each period of 
the modulation. The hrst sweet spot of this band, at low T, corresponds to solutions that 
complete one positive and no negative phase slip per cycle; the second sweet spot, further 
up, is comprised of solutions that complete two positive and one negative phase slips per 
cycle, etc. The second region on the right corresponds to solutions that undergo two extra 
positive phase slips, and so on as vq increases. All these regions have a similar structure 
as the modulation period T increases. They all correspond to the resonance tongues in 
Fig. [island are separated by transition zones with solutions that have a non-integer winding 
number. These transition zones narrow as T increases and solutions within them can have 
periods that are a multiple of the modulation period, or not be periodic at all. Solutions with 
negative winding number are found in analogous regions obtained by reflection in ro = 0. 

Figure [15] shows the winding number between PO and the 1:1 resonance tongue as com¬ 
puted from time simulations averaged over 5000 modulation periods T = 25. The hgure 
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FIG. 15. The winding number N as a function of tq across the transition zone between PO and 
the 1:1 resonance tongue for T = 15 and a = 2. 


shows that the winding number increases monotonically and smoothly within this transition 
zone, as expected on the basis of theoretical considerations [^. However, modihcations 
of the nonautonomous Adler equation, such as the inclusion of an inertial term or a more 


genera. 

zones 


time dependence, can generate subharmonic resonances that populate the transition 


16| . Subharmonic resonances have also been observed to produce a devil’s staircase 


type structure in the related problem of spatially localized states in the periodically forced 
Swift-Hohenberg equation 241. 


C. Asymptotic formation of sweet spots 

We can extend the above predictions by analyzing a limit in which the trajectory barely 
exits the phase-locking region — 1 < r < 1 but assuming the modulation period is slow 
enough that phase slips still take place. Explicitly, we take r(t) = -|- (1 -b e^p) sine^cut, 

where e^/i represents a small offset of the average value of r{t) from tq = 0. We introduce 
the slow time scales r = et and <h = e^ut and employ an asymptotic expansion of the form 
9 = 9o + eOi + 6^02 +- 

At leading order, the Adler equation ([1]) gives sin 9o = sin $ for which we choose the 
stable phase locked solution 00 = *^-!- 27m that has no r dependence. The alternate choice. 
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6*0 = TT — $ + 2 ?™, produces unstable periodic orbits or unstable winding trajectories. At 
order e, we obtain the equation 9^6*o = — 6 *i cos 6 *o. When 6 *o 7 ^ 7 r /2 + vm, 6*1 = 0 in order to 
satisfy the condition that 6*0 be independent of r. At order e^, we obtain 

6*2 cos 6*0 = /i + p sin <h + ^ 6 *^ sin 6q — drOi — a;9#6*o, (34) 

leading to the second order correction 

6*2 = (p — 1 ^) sec <h + ptan $ (35) 


provided that 6*0 7 ^ 7r/2 + Tin. 

To examine the dynamics near 6*0 = 7 r /2 + nii where the system is transitioning between 
phase-locked dynamics and winding, we take the slow time to be $ = 7 r /2 -|- e0. Equation 
(fT|) then becomes 

eud(j)9 = e^p + (1 + e^p) cos ecj) — sin 6. (36) 

The leading order and order e equations are identical to the general case above while 6*1 is 
determined from the order equation, 


6*2 cos 6*0 = p -6 p — sin 6*0 — ujd(j,9i, (37) 

which differs from Eq. (l34ll . Since 9o = 7 r /2 the Riccati transformation 9i = —2u:d^'ipl'ip 
transforms this equation into the Weber equation 

+ (38) 


We now use a matching procedure to connect the relevant solution of this equation to the 
case when $ 7 ^ 7 r /2 -|- nn. Noting that 6*i(<h) —)■ 0 as —)■ 7 r /2 -|- rni, we choose our solution 
such that 6*1 —)■ 0 as 0 ^ — 00 . This matching condition is satisfied by the parabolic cylinder 
function -0 = Dy(s), where 


p -6 p 1 _ (f) 

2u 2 ’ ^ 


(39) 


Each zero s = sq of V’ = Du(s) corresponds to one phase slip. Care must be taken in 
interpreting these results since the zeros of t/t correspond to divergences of 6*1 and thus a 
breakdown of the asymptotic series used to obtain Eq. (1381) . The above calculation holds 
between the asymptotic breakdowns where 6*1 (r) diverges, so a complete trajectory can be 
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constructed by “gluing” solutions across each individual phase slip. Thus ip can be used to 
describe a series of phase slips via this gluing process. 

The number of zeros of pj corresponds to the number of phase slips and thus determines 
which solution branch the system will follow upon re-entering the phase-locked region 7i/2 < 
<h < 37r/2. In particular, [n+] phase slips are undergone when 

M “ I < M + i (40) 


More generally, we can express the number of positive (negative) phase slips that occur 
near the boundaries of the phase-locked region in terms of the parameters of the problem as 




{a±ro-l)T 

An 


(a ± ro — 1) > 0 

(a ± ro — 1) < 0, 


(41) 


where the square bracket indicates rounding to the nearest integer. The predictions of 
this theory match well with time simulations for a = 1.005, as seen in Fig. [161 The 
simulations employed a fourth order Runge-Kutta scheme for 12 periods of the modula¬ 
tion using the initial condition 0(0) = sin“^ tq. The winding number was computed from 
27iN = (0(12T) — 6{2T)) /lO as a function of the parameters vq and T. Time simulations 
were used in place of numerical continuation because the extremely long time scales make 
continuation a computationally challenging task. Owing to symmetry the PO region is al¬ 
ways centered on tq = 0, and states with negative winding number are found in regions 
obtained by reflection in ro = 0. 

The hgure reveals the formation of sweet spots in this limit whenever a > 1. When 
a < 1, there are two distinct sets of resonance bands - one set formed by regions with 
a hxed number of positive phase slips n+, and the other by regions with a hxed number 
of negative phase slips n_. At a = 1 the two sets of resonance bands both asymptote 
to ro = 0 as T oo (Fig. ITTIf a)). The sweet spots and pinched zones emerge through 
the intersections of these resonance bands that take place once a > 1 (Fig. ITW b.c)). In 
particular, the pinched zone separating the n and n -|- 1 sweet spots in the PO region is 
located at (a — l)T/47r = n -|- 1/2 and marks the transition from n to n -|- 1 positive and 
negative phase slips within a modulation cycle. 
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FIG. 16. (Color online) Average winding number per period T of the frequency parameter shown 
in the [tq^T) plane for a = 1.005. Colors represent results from numerical simulation: no net 
phase slips occur over the course of a modulation period in the dark region {PO) to the left; the 
alternating yellow/orange regions to the right indicate 1,2,3,... net phase slips as vq increases. 
The red (negative slope) and blue (positive slope) lines mark the transitions between the regions of 
constant and n_ as predicted by the asymptotic theory (Eq. (I^T]) h Trajectories with negative 
winding number are located in regions obtained by reflection in rg = 0. 

V. ADIABATIC THEORY 

We now consider a more general time-dependence for the parameter r, but assume it 
varies slowly enough that we can treat the dynamics quasi-statically: r = r{2nt/T) with 
T !:§> 1. In this adiabatic limit, two distinct types of dynamics arise: slow dynamics that 
track the steady state phase-locked solution when — 1 < ri^nt/T) < 1 and a fast phase 
rotation with an adiabatically varying parameter when |r(27rt/T)| > 1. No matter how 
low the frequency is, there is always an intermediate regime around the transition from a 
phase-locked state to rotation where the phase rotation is slow enough that it occurs on the 
same scale as the parameter drift. We apply WKB theory to capture the dynamics in each 
of the two regions separately and provide a condition for matching the solution across the 
transitions at r ^ ±1. 

We start with Eq. fl32|) but assume that r = r{ojt) with cu <C 1. We do not need to specify 
the form of r{ut). We transform this equation into a standard form for WKB theory by 




27 



FIG. 17. (Color online) Transitions between the regions of constant n+ (red, negative slope) and 
n_ (blue, positive slope) in the {ro,T) plane as predicted by the asymptotic theory (Eq. (HD) ') for 
a = 1.0000, 1.0025, 1.0050, respectively. A sweet spot and pinching structure begins to emerge as 
a increases. 


recasting it in terms of the slow time 0 = ut\ 


'< I 

^ H- 


— 1 car' 

- "a -^ "o-^ 

4 2r 


uP'r" 

2r 


3u\r'f 


z = 0. 


(42) 


The system transitions from a phase-locked state to winding near — 1 ~ 0{u}), and we 
can use the standard WKB ansatz 2 ; = A2;wkb = Aexp(zS'/a;) -f c.c., where A is an arbitrary 
complex constant determined from initial conditions and/or matching procedures when we 
are away from these points. We suppose that S = So + uSi + ... and match orders to solve 
for each Si. 

Making the WKB substitution generates, at leading order. 


q/2 _ 

^0 — 


— 1 


(43) 


The leading order WKB solution, in terms of the original time scale, is z = A exp (±| / \/r^ — Idtj. 
The equation at next order is, after simplihcation. 


yielding 


2 ^'+ tS'^ + ^ = 0 , 


Si = - { i log \/r^ — 1 =F tan ^ 


\/r2 — 1 


(44) 


(45) 


depending on the choice of root for So- 

Including this correction, the solution becomes 

A 


1)V4 


exp ±1 ( / Vr‘^ — Idt — tan ^ 


Vr 


2; = 


(46) 
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When r < 1, we find it convenient to rewrite the expression for Si as 


Si = ^ flog Vl - r2 ± log 


1 +vT 


(47) 


and the solution now takes the form 


A 


z = 


1 +vT 


T1/2 


exp / Vl — r‘^dt. 


(48) 


(1 -r2)V4 

Near a transition point r = 1 (we take it to be at t = 0), we suppose that r ^ 1 + at, 
where a = r(0) = a;r'(0) is a constant. To leading order, the equation becomes 


a 


z + -{t+l)z = 0, 


(49) 


which has solutions in terms of the Airy functions Ai(s) and Bi(s), where s = — (^+1)- 

We will further assume a > 0 so that the transition occurs as the system leaves the 
phase-locked region and enters the winding region, and remind the reader that 

1 


9 2z 

tan - =-h - 

2 r z r 


1 - 


(50) 


We consider the solution within the phase-locked region that follows, for f < 0, the stable 
steady-state solution branch 9 = sin“^ r, corresponding to taking the negative root of S'q. 
Thus, equation fl45]l reduces to 


A 


z = 


pi 


(1_ ^2)1/4 


1 + ^r^ 


exp i / Vl — r'^dt, 


(51) 


where Api depends on the choice of initial condition. In terms of 9, expression IHTl reads 


tan - = 
2 


9 1 - 


(52) 


r(l - r2) 

In order to match solutions across the transition region, we must take the t —)■ 0 (equivalently 
r —)• 1 limit of this solution and match it to the t —)• —cxo (equivalently s —>■ cxo) limit of 
the Airy function solution of Equation (09]). This procedure selects the Airy function Ai 
with amplitude proportional to Api. On the other side, the winding solution coming from 
equation (061) can be matched to the Airy solution when written in the form 

Aw 


r 


2 _ 1)1/4 CO® 2 


\/r2 — \dt — tan ^ 


= 


■\/r2 — 1 


(53) 




















29 



FIG. 18. (Color online) The phase 6{t) mod(27r) near the transition from a phase-locked state to 
winding from a time simulation with T = 27r x 103, a = 2, and ro = 0 (black solid line). The 
simulation represents the evolution of 9 in the time window [—100; 160] of a converged periodic 
orbit for which t = 0 corresponds to r = 1. The dashed lines are computed using the adiabatic 
predictions (I52p . (I54p without the “subdominant” term proportional to r while the dotted lines 
take it into account. Predictions in the phase-locked (winding) regime are shown in blue (green) 
for t < 0 {t > 0). 


where = ^^(Api). The matching is achieved by comparing the r ^ 1 (t —)■ 0) limit of 
expression fl5^ to the t ^ oo {s ^ —oo) limit of the Airy function obtained in the matching 
procedure with the phase-locked solutions. Expression fl53|) yields: 


9 1 

tan - = - 
2 r 


1 -|- -v/r^ — 1 tan ■ 


\/r^ — Idt — tan ^ 


1 - 


r(r2 — 1) 


(54) 


Vr^ — 1 , 

Figure [18] shows a comparison of the WKB solution in terms of 6 with a periodic orbit 
obtained through simulation with r{t) = 2sin(10“^f -|- vr/6). 

The results obtained from the WKB approximation in the limit of a slowly-varying fre¬ 
quency parameter can be generalized using a theorem that places bounds on the number 
of zeros of solutions to linear second order differential equations. Given an equation of the 
form 

z -I- q{t)z = 0 (55) 


with q{t) > 0 in and bounded, such that q{t) = o {q^^‘^{t)) as t —>■ cx), it can be shown 25] 


















































30 


that the number of zeros [n] between 0 < t < T for a given solution z{t) 7 ^ 0 is bounded by 


7r[n] — 

It follows that when q <^1 


\/q{t)dt 


< TT + 




16g5/2 4g3/2 


dt. 


7r[n] ~ / \/q{t)dt as T —)■ 00 , 


(56) 


(67) 


thereby reproducing the quasi-static prediction from WKB theory. In the case of the Adler 
equation, the corresponding frequency parameter is given by 


— 1 r r 3r^ 
^ 4 2r 2r 4r^ 


(58) 


The conditions on r for the applicability of the bound within the time interval of interest 
are that \q{t)\ > 0, r G and is bounded. We can make some further approximations in 
the limit that r = r{ut) is slowly varying, i.e., a; 1 , and the hrst condition reduces to 
|r| -|- 0{uj) > 1. In this adiabatic limit, the integral in the bound becomes 

1 


\/q(f)dt = 


\/r^ 


-dt — tan 


-1 


- 1 


+ 0{uj). 


( 69 ) 


The bound on the number of zeros of the Hill equation translates into a bound on the 
number of phase slips incurred by a solution to the Adler equation over a given time interval 
where q{t) > 0, i.e., when r{t) is outside of the phase-locking region. We dehne n± by the 
integral 

n± = — [ \/q(t)dt (60) 

^ Jt± 

over the time interval T± spent with q{t) > 0 and r{t) > 1 (r(f) < —1) for ? 7 ,+ (n_). The 
bound described above restricts the number of phase slips over 7± to either rounding up 
or down ([n±J or [n±]) to order 0{uj). This is a generalization of the WKB solution in 
the sense that the bound applies even when the slowly-varying assumption does not hold. 
Some care must be taken when applying this bound as g —>■ cx) as r —>■ 0. The bound must 
be applied to positive and negative phase slips separately in order to place a bound on the 
winding number of a particular trajectory. 

The WKB approximation can be used to predict the partitioning of the parameter space 
by winding number (see Fig. fMl) by computing the net winding number N = [n+] -f [n_], 
where 

T ^- 

n± = ±—J v(^o ± asin0)2 - ld0. 


(61) 
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(a) WKB 



(b) bound 



FIG. 19. (Color online) Average winding number per period T of the frequency parameter shown in 
the (ro,T) plane for a = 2. Colors represent results from numerical simulation: no net phase slips 
occur over the course of a modulation period in the dark region to the left; the alternating lighter 
yellow/darker orange regions to the right indicate 1,2,3,... net phase slips as tq increases. The 
red/blue (negative/positive slope) lines represent predictions of adiabatic theory. The left panel 
shows the prediction based on the WKB approximation (Eq. (1611) 1 while the right panel shows the 
prediction based on the bound in Eq. (|57)) . 

and tq ± asin0± = 1. The first correction from WKB theory cancels becanse the system 
always enters and exits the phase-locked region at the same valne of r. Replacing the 
expression in the sqnare root with q{t) provides a way to estimate the winding nnmber from 
the bonnd. Figure [19] shows a comparison of the resulting prediction with the numerical 
results in Fig. [TH We see that the adiabatic theory agrees well with the numerical results 
far beyond the low frequency limit for which it was constructed, a conclusion supported by 
the generalization (jSSj). 


VI. DISCUSSION 


In this paper, we have investigated the dynamics of two coupled oscillators when the 
frequency difference is modulated in time. The same equation describes a multitude of other 
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systems, ranging from Josephson junctions to systems of large numbers of coupled oscillators 


as detailed in Sec. [H Specifically, we studied here the Adler equation [l| with a sinusoidally 
varying frequency parameter. The frequency modulation introduces two new parameters 
into the problem, in addition to the mean frequency difference tq: the amplitude a and the 
period T of the modulation. While the autonomous Adler equation leads to phase locking 
for — 1 < ro < 1 and persistent drift for |ro| > 1, we have unveiled much richer dynamics 
that take place when frequency modulation is activated: the phase-locked solutions turn 
into periodic orbits and the phase difference 6 between the oscillators becomes a periodic 
function of time. The region PO of the existence of these periodic orbits is centered around 
ro = 0 and exhibits a succession of sweet spots as a or T increases, interspersed with pinched 
zones where the width of the PO region vanishes. The width of these sweet spots decreases 
with increasing a and T. On either side of PO are regions within which the solution grows 
or decays by one, two, etc. phase slips per modulation cycle. These regions have the same 
basic structure as the PO region and are separated by exponentially thin transition zones 
where the number of phase slips fluctuates from cycle to cycle. This intricate behavior is a 
consequence of a sequence of resonances between the time needed for a phase slip and the 
period of the modulation, and can be described, in an appropriate regime, in terms of an 
interaction between n:l and —n:l resonance tongues. 

Canard orbits form an essential part of this picture 2q. These are present in the vicinity 
of the boundaries of the PO region and consist of trajectories that drift along a branch 
of stable equilibria for part of the cycle; after reaching a fold at which the equilibria lose 
stability the trajectory drifts for a time along the branch of unstable equilibria, instead of 
detaching, before an abrupt jump back to a stable equilibrium. Equation fl38|l describes the 
emergence of such trajectories for low frequency modulation with mean near ro = 0 and 
amplitude slightly larger than a = 1; Fig. 1201 shows several examples of the predicted canard 
solutions, for comparison with the “larger” periodic canard orbits computed numerically in 
Sec. HTTEI 

We mention that similar behavior has been observed in the partial differential equation 
description of the dynamics of spatially localized states 23|]. In this work, the quadratic- 
cubic Swift-Hohenberg equation (SHE23) is forced in a time-periodic manner and a similar 
partitioning of parameter space is observed (Fig. [2T]) . The reason for this similarity can 
be traced to the nature of the motion, under parametric forcing, of fronts connecting a 
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FIG. 20. (Color online) Canard behavior near r = 1 in the limit T 1 as predicted by Eq. (j38p 
for V = —10“^ (red, inner), —10“® (blue, middle), —10“^^ (green, outer) and 0 (black). In terms of 
the parameters of the original problem v = + a — 1) — the horizontal and vertical scales 

are r — 1 ~ 1/T and 9 — 7r/2 ~ l/\/r. The stable (solid purple) and unstable (dashed brown) 
stationary solutions to the autonomous problem are shown for reference. 


spatially periodic state of SHE23 to the trivial, homogeneous state: the front motion is 
analogous to repeated phase slips, with each “phase slip” corresponding to a nucleation or 
annihilation event that adds or subtracts one wavelength of the pattern at either end of 
the localized structure. However, the resulting partitioning of the parameter space is not 
symmetric owing to a lack of symmetry between positive and negative “phase slips”. An 
adiabatic theory of the type described here works equally well in SHE23 and its predictions 
are in excellent agreement with the results of numerical simulations 2^. Indeed SHE23 
also displays canards associated with the transitions between different states (lightest gray 


regions in Fig. 


m- 


The work presented here has a direct application to Josephson junctions driven by an 
AC current. In the overdamped limit such junctions are modeled by Eq. ([T]) 1^, with r 
representing the external drive and 9 the phase difference of the Ginzburg-Landau order 
parameter across the gap. The so-called supercurrent across the gap is proportional to sin 9 
while the voltage produced corresponds to the derivative 9. In this context, phase-locking 
and phase-slips are closely related to the existence of Shapiro steps ^ for a single Josephson 










34 



FIG. 21. (Color online) Average winding number per period T shown in the {rQ,T) plane for the 
Swift-Hohenberg equation when 6 = 1.8 3, 24 1. In the alternating lighter yellow and darker 
orange regions on the right, a localized state grows by the net addition of 2,4, 6,... wavelengths 
of the pattern within each forcing cycle, while in the alternating light and dark blue regions on the 
left it shrinks by 2,4,6,... wavelengths within each cycle; the dark region labeled PO corresponds 
to localized states that pulsate but maintain constant average length. The gray areas represent 
transition zones where the average winding number per period is not an integer. In this system, 
such zones are characterized by a devil’s staircase type structure. 


junction. Related dynamics arise in arrays of Josephson junctions that are globally coupled 
via an LRC circuit 2l|]. These systems provide a physical realization of the phase-coupled 
oscillator models mentioned in the introduction. 


In fact, weakly coupled systems can often be decomposed into two parts with part A 
obeying dynamics that are largely insensitive to the dynamics of part B. In these circum¬ 
stances it often suffices to consider system B on its own but with prescribed time-dependence 
arising from the coupling to A. This is the case, for example, in globally coupled phase os¬ 
cillator systems, in which each oscillator responds to the global dynamics of the system but 
the global dynamics are insensitive to the details of the dynamics of an individual oscillator. 
These systems, for reasons explained in the introduction, have properties closely related to 
the nonautomous Adler equation studied here. For these reasons we anticipate applications 
of the techniques developed here to studies of synchronization in oscillator networks. 
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